T • H • E 


OHIO 

SIXIE 

UNIVERSITY 


Z. t-S 




/ As -c/2^ 


7/v/' 


Analysis of Lossy Composite Terminating Structures 


R. Andre, A. Dominek, J. Munk and N. Wang 


The Ohio State University 

ElectroScience Laboratory 

Department of Electrical Engineering 
Columbus, Ohio 43212 


Technical Report 723224-3 
Grant No. NAG3-1000 
September 1991 


National Aeronautics and Space Administration 
Lewis Research Center 
21000 Brookpark Rd. 

Cleveland, OH 44135 


(NASA-CR- 189901 ) ANALYSIS OF LOSSY 
COMPOSITE TERMINATING STRUCTURES (Ohio 
State Univ.) 69 p CSCL 20N 


G3/32 


N92-1921 7 


Unci as 
0071411 


A. Approved for public release; Distribution is unlimited 


NOTICES 


When Government drawings, specifications, or other data are 
used for any purpose other than in connection with a definitely 
related Government procurement operation, the United States 
Government thereby incurs no responsibility nor any obligation 
whatsoever, and the fact that the Government may have formulated, 
furnished, or in any way supplied the said drawings, specifications, 
or other data, is not to be regarded by implication or otherwise as 
in any manner licensing the holder or any other person or corporation, 
or conveying any rights or permission to manufacture, use, or sell 
any patented invention that may in any way be related thereto. 


502T2-101 


REPORT DOCUMENTATION 
PAGE 

1. REPORT NO. 

2. 

3. Recipient's Accession No. 

4. Title and Subtitle 


/ 


6: Report Date 

September 1991 

Analysis of Lossy Composite Terminating Structures 



e. 

7. Author(s) 

R. Andre, A. Dominek, J. Munk and N. Wang 

8. Performing Org. Rept. No. 

723224-3 

9. Performing Organisation Name and Address 

The Ohio State University 



10. Project/Task/Work Unit No. 

ElectroScience Laboratory 
1320 Kinnear Road 
Columbus, OH 43212 




11. Contract(C) or Grant(G) No. 
(C) 

(G) NAG3-1000 

12. Sponsoring Organization Name and Address 

National Aeronautics and Space Administration 



18. Report Type/Period Covered 

Technical Report 

Lewis Research Center 

21000 Brookpark Rd., Cleveland, 

Ohio 44135 



14. 

15. Supplementary Notes 

IB. Abstract (Limit: 200 words) 






A finite element solution and computer code for the electromagnetic scattering of inhomogeneous 
penetrable bodies is presented. The application for the code is for the analysis and design of leading 
and trailing edge terminations when conducting and non-conducting materials are used. Examples of 
simple triangular shaped terminations are also presented. 

17. Document Analysis a. Descriptors 
MATERIALS 
SCATTERING 
MOMENT METHOD 
b. Identiflers/Open-Ended Terms 


RCS 

2-D 



c. COSATI Field/Group 






18. Availability Statement 

A. Approved for public release; 


19. Security Class (This Report) 

Unclassified 

21. No. of Pages 

69 

Distribution is unlimited. 


20. Security Class (This Page) 

Unclassified 

22. Price 


(See ANSI-Z39.18) 


See Instructions on Reverse 


OPTIONAL FORM 272 (4-77) 
Department of Commerce 


i 
























Contents 


List of Figures v 

1 Introduction 1 

2 Theoretical Analysis 3 

I Interior Region 5 

• II Exterior Region 10 

II. 1 Absorbing Boundary Method 10 

II. 2 Hybrid FEM/BEM 14 

III Matrix Inversion 23 

IV Scattered Far Fields 26 

V Geometry Specification : 27 

3 Examples 29 

I Validation 29 

II Terminations 33 

4 Conclusions 44 

iii 

PRECEDING PAGE BLANK NOT FILMED 



Appendix 


A BEM/FEM Numerical Implementation 45 

I Introduction 45 

II Geometry 45 

III Component Tn^* ei 51 

IV Component m^. e i 55 

V Component v e ,i 59 

VI Far Field 60 


IV 



List of Figures 

2.1 Illustration of conceptual edge termination 4 

2.2 Typical first and second order triangular elements 8 

2.3 Matrix distribution of non-zero elements before renumbering. . 24 

2.4 Matrix distribution of non-zero elements after renumbering. . 25 

3.1 Test geometries: (a) Concentric cylinders, (b) Flat strip. ... 31 

3.2 Mesh used for the concentric cylinders 32 

3.3 Mesh used for the flat strip. The centerline of the dashed 

region is the location of the strip 33 

3.4 TM bistatic echo width for the dielectric concentric cylinders. 

FEM/BEM - solid, Eigenfunction - dashed 34 

3.5 TE bistatic echo width for the dielectric concentric cylinders. 

FEM/BEM - solid, Eigenfunction - dashed 35 

3.6 TM bistatic echo width for the flat strip. FEM/BEM - solid, 

Eigenfunction - dashed 36 

3.7 TE bistatic echo width for the flat strip. FEM/BEM - solid, 

Eigenfunction - dashed '. 37 

3.8 Test terminations, (a) Terminations Tl and T2, (b) Termina- 
tion T3 38 

3.9 Echo widths for termination Tl. TM - solid, TE - dashed. . . 39 

3.10 Echo widths for termination T2. TM - solid, TE - dashed. . . 40 

3.11 Echo widths for termination T3. TM - solid, TE - dashed. . . 41 


v 



3.12 The internal TE scattered field for termination T3 with illu- 
mination from 25°. 42 

A.l A five sided boundary. The numbering of the nodes and 

boundary elements is shown 46 

A. 2 Linear generating functions on the boundary elements 47 

A.3 Elemental linear generating functions on a single boundary 

element 48 

A.4 Quadratic generating functions on the boundary elements. . . 48 

A. 5 Elemental quadratic generating functions on a single boundary 

element 49 

A. 6 The division of the unit square into upper and lower triangles. 58 


vi 



Chapter 1 
Introduction 


Leading and trailing edges on aerodynamic surfaces often have special edge 
terminations due to design constraints. One constraint is their impact upon 
the electromagnetic scattering performance for aerodynamic surfaces. The 
edges are important in determining the scattered field from these surfaces. 
Their impact is seen in several different ways. For an electrically thin, per- 
fectly conducting surface, near grazing incidence, the leading edge is primar- 
ily sensitive to an electric field polarized parallel to the edge while the trailing 
edge is primarily sensitive to an electric field polarized perpendicular to the 
trailing edge. 

The previously mentioned scattering behavior changes as the surface be- 
comes electrically thicker. The leading edge becomes a specular scatterer and 
is sensitive to both polarizations while the trailing edge will allow energy to 
“creep” around the edge. Other changes occur when the surface is no longer 
perfectly conducting due to fabrication considerations. Now the incident field 
can penetrate the structure and scatter with greater complexity. 

The character of the scattered field depends upon several parameters. The 
most important parameter is the shape of the structure. The scattered field 
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is focused in directions where Snell’s law is satisfied. Backscattered fields are 
the strongest where the incident field direction is normal to a flat surface. A 
curved surface results in lower backscattered fields but then provides broader 
angular scattering. Hence, it is important to properly shape both the exterior 
and interior of penetrable surfaces. 

Another parameter involves the various materials used in a penetrable 
surface. Abrupt changes in material with different electrical properties can 
also contribute to the scattered field. Thus, it is then desirable to taper 
the electrical material properties when possible. Appropriately chosen loss 
characteristics can also help in reducing the internal visibility of a structure 
when sufficient tapering is not possible. 

Hence, the design of terminating structures involves the judicious choice 
of both structural shape and material components to achieve the desired 
electromagnetic scattering performance. A computer program is under de- 
velopment to assist in the design of two-dimensional edge terminations. This 
code is based upon the rigorous solution of Maxwell’s equations. 

The following chapters discuss the general formulation used in the code, 
its present status and a few examples to demonstrate the validity of the 
code. A future report will present the final form of the code, its operation and 
comparisons with measurements to illustrate particular design considerations 
for edge terminations. 
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Chapter 2 

Theoretical Analysis 


The analysis chosen for the basis of the computer program to calculate the 
scattered field from an edge termination was based upon its potential ge- 
ometry and its material properties. Figure 2.1 illustrates a conceptual ter- 
mination which has several homogeneous regions. Each region can be de- 
scribed geometrically with unique constitutive parameters. There are two 
basic approaches to solve such problems. The traditional approach involves 
the formulation of an integral equation with unknown equivalent volumetric 
currents occupying the space of the scatterer [1]. This approach is called a 
method of moments solution (MM) where the integral equation is discretized 
into a system of equations. 

The alternate approach involves the direct solution of the differential form 
of Maxwell’s equations [2, 3]. This approach is called a finite element solution. 
Traditional finite element solutions (FEM) have been employed in problems 
with bounded (finite) regions but their presence is becoming more common 
in problems with unbounded regions such as in electromagnetic scattering 
problems where the region extends to infinity. 
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Figure 2.1: Illustration of conceptual edge termination. 
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There are several differences between the two basic approaches. The most 
noticeable difference involves the final matrix for inversion to solve for the 
unknowns. The MM solution yields a dense matrix while the FEM solution 
yields a sparse matrix. This difference becomes significant when the problem 
becomes electrically large since the number of unknowns will also increase. It 
is for this reason that the FEM approach was taken to develop the necessary 
scattering code. 

The following sections briefly describe the FEM formulation and its im- 
plementation. The work is centered at determining the local and far fields 
from and inhomogeneous scatterer illuminated by a plane wave. The analysis 
is divided into two parts, i.e., the bounded internal region which consists of 
the scatterer and the unbounded external region which includes all of space 
excluding the scatterer. 

I Interior Region 

There are two independent scattered field solutions to Maxwell’s equations 
for two dimensional geometries. Defining the geometry cross-section in the 
x-y plane, the solution is determined by the polarization of the incident field. 
The solution is called TM (transverse magnetic) if the incident electric field 
is z polarized or TE (traverse electric) if the incident magnetic field is z 
polarized. The scalar wave equation to be solved for the total field where the 
relative permittivity (e r ) and permeability (fi r ) vary with position is 

V ‘ (“ V< ^) + \ k2<i> = ° (2.1) 
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where 



TM polarization. 


TE polarization. 


( 2 . 2 ) 


(2.3) 


with A: being the propagation constant in the medium which equals ^ Jt7\L r • 
The total field is the sum of the incident and scattered fields denoted with 
superscripts i and a , respectively. 


The total field solution can be obtained through using the technique of 
weighted residuals. This technique involves the integration of the differen- 
tial equation over the bounded region S times a testing function, wj. The 
resulting expression from this operation is 


I LI dxiy = L 


(2.4) 


where we have made use of the vector identity 


WjV • = V • 


(2.5) 


The line integral is a result of applying the two dimensional divergence the- 
orem 'to simplify Equation (2.4). This boundary term associated with 9T, is 
dependent upon certain boundary conditions discussed later. 

The next step is to expand (f) in some basis and for simplicity, the same 
expansion as the testing function, Wj , is used for both. The expansion used 
are the common triangular elements [2] although there are several others. 
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The basis is a polynomial approximation represented by subsectional ele- 
ments which span the bounded geometry region. The expanded region is 
represented with individual elements as shown below 


= '52M x ,y) (2-6) 

e=l 

with each element, <f> e (x,y)\ being represented as 

N„ 

(2.7) 

1=1 

where tn, are polynomial shape functions, N e are the number of elements and 
<j>i represents the field at the tth node of the eth element. The constant N„ 
equals 3 or 6 depending upon if an element consists of three or six nodes as 
shown in Figure 2.2. A three nodal element corresponds to a linear approx- 
imation to the fields while a six nodal element allows a quadratic variation 
in the fields. 

It is convenient to define the test functions such that w, = 0 at all nodes 
except one, where it has unity value. This is accomplished by letting 

Wi = Ri(n,ti)Rj(n,£ 2 )Rk(n;t3), i+j + k = n (2.8) 

where 

1 m — 1 

Rm(n,() = -]!("(■ *)> m > 0 

m - *= 0 

J2o(»,0 = 1. (2.9) 

The interpolation functions wj may be expressed in terms of simplex 
coordinates which are independent of the global x — y coordinates to simplify 
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the necessary calculations [2]. The relationship between the cartesian and 
simplex coordinates is given by 

6 1 j [ *2i/3 - *3j/2 y 2-1/3 *3 - *2 | [ 1 

6 ~ Ya X3Vl ~ XiV3 V3 ~ yi Xl ~ X3 x (2.10) 

. 6 J x lV 2 - X 2 V\ 2/1 - 1/2 *2 - *1 J y . 

where 

i 1 zi yi ' 

A = -det 1 * 2 1/2 ( 2 -ll) 

1 *3 1/3 . 

represents the area and y\,x 2 etc... are the vertices of the element under 
consideration. 

In matrix form, Equation (2.4), can be written as 

Se = v (2.12) 

where 

s„ = JJ sZ (v„,.v w , — k 2 W{Wj ) dxdy (2-13) 

and 

Vi = — f Wih ■ V<j>dl (2.14) 

v Jd r 

Efficient evaluation of the terms in Equation (2.13) is provided in [2]. 

The final boundary consideration applies to conducting boundaries. The 
nodes along a conducting boundary are known for the TM case since E l z = 0 
and are excluded from the matrix equation as unknowns. For the TE case, 
since the normal gradient is identically equal to zero along the conductor, 
the nodal points along the boundary are treated as unknowns. 
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II Exterior Region 


The FEM approach allows a very convenient boundary condition implemen- 
tation for bounded regions. However in scattering problems, the region is 
unbounded (boundary at infinity) and a far field boundary condition has to 
be employed. There are two basic approaches to accomplish this task. One 
approach is to construct a boundary near the scatterer which ideally absorbs 
all outward propagating waves and creates no reflected waves. This approach 
is called the absorbing boundary condition (ABC). 

The second approach couples the fields calculated by a FEM calculation 
(bounded region) to the exterior through a integral expression provided by 
the boundary element method (BEM). This approach is called the hybrid 
FEM-BEM approach and is theoretically more accurate since no scattered 
field approximations are made at the coupling boundary as is done for the 
ABC method. 

II. 1 Absorbing Boundary Method 

The goal in this technique is to make a bounded region from an unbounded 
region in order to directly use the FEM. This can be accomplished only ap- 
proximately since the Sommerfeld far field radiation condition is only asymp- 
totically satisfied at some boundary which encloses the scatterer. The Som- 
merfeld far field radiation condition 

lim -- = —jk<f> a (2.15) 

p—oo op 

is valid for 0r — > oo. However, for numerical analysis it is necessary to 
apply a boundary condition valid for 9T located at a finite distance from the 
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scatterer. From a practical standpoint, it should be as close as possible but 
far enough away for the following approximations to be valid. 


Expanding the scattered field asymptotically 


Mp,o) = 


o-jkp 


Vp 


|<£o 


(*) + 


MO) , M^) , MO) 


+ 


+ 


+ 


(2.16) 


Imposing the Sommerfeld far field radiation condition on c?r results in only 
the first term of the asymptotic expansion being satisfied. Bayliss and Turkel 
[5] presented a higher order radiation condition 


% = 


which satisfies the first four terms in Equation (2.16), where 


(2.17) 


a(p) = 


~ T P ~ g & 




"1 


(2.18) 


and 


P(p) = 


1 


2ki? 

j- 
kp 


(2.19) 


Substituting Equation (2.17) into Equation (2.4) with <p s = <f> — <j>' and 
integrating-by-parts yields 


J J — (Vwj • V<f> — K 2 wj<j)) dxdy 

-L( aWi +-^ s $ pM 

= L Wi {w - _ p w) pie 


( 2 . 20 ) 


where dr is circular and centered at the origin. In Equation (2.20), <t>' is 
assumed known while <f> is the unknown variable. 
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In matrix form, Equation (2.20), can be written as 


Se = v. 


( 2 . 21 ) 


The matrix elements are given by 

s* = 5 ,; + s'j 


( 2 . 22 ) 


where 


S-j = J J — (Vw, • Vwj — k 2 W{Wj ) dxdy 


(2.23) 




aw{Wj 


dwi dwj ^ 

He He ) 


pd9 


and 



-00-0 


d 2 4>' 
dd 2 


t) 


pdd. 


(2.24) 


(2.25) 


The terms expressed in Equation (2.23) have been previously computed, thus 
only the second integral need be considered. Due to limitations on the num- 
ber of triangles one can expect to accommodate, approximating the curved 
boundary £r by a straight line introduces discretization error. If, however, we 
do not restrict our trial function polynomials to exist only within the triangle 
but allow them to “spill” over onto the curved boundary we can eliminate 
this error [6]. This is valid since the polynomials exist over all space. The 
node values however, are still given at points on the triangle. In polar form 
the global coordinates are 


\ti 

(2 

. (3 . 


' ® 2 y.3 - *3j/2 V2 - y3 *3 - *2 " 


1 ' 

*3j/l - V3 - Hi *1 ~ *3 


p cos 9 

xii /2 - * 2 yi yi -V 2 *2 - *1 . 


p sin 6 


(2.26) 
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The test functions used in evaluating the integrals are the same functions 
previously given in Equation (2.8). Additionally, the partial derivatives in 
Equation (2.24) are expanded in terms of the area coordinates which yields 

di Vi _ JLv dwi d(i 

~d0 

Equation (2.24) is valid for a circular absorbing boundary centered at the 
origin of our coordinate system, therefore p is constant and Equation (2.24) 
can be solved in closed form. 

Finally, the voltage coefficients v, given in Equation (2.25) are considered. 
The incident plane wave with a e~* wt time convention is used and is defined 
as 



<t>\p,6) = «£ 0 e^ cos( *- 7) 


(2.28) 


where <j) u is the magnitude of the plane wave and 7 is the angle of incidence. 
Taking the derivatives in Equation (2.25) results in 


d(f>' 

dp 


= jkcos(6 — 7)<^'(p,0) 


(2.29) 


and 


^7 = {-jkpcos(6 - 7) - k 2 p 2 sin 2 (0 - 7)} (j>'(p,e ). (2.30) 

Inserting Equations (2.29) and (2.30) into Equation (2.25) and numerically 
integrating yields the required coefficients, v,. 
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II. 2 Hybrid FEM/BEM 


The efficiency of the FEM technique for unbounded problems can be greatly 
increased if the required meshed region can be reduced. The previously 
presented ABC technique allowed a convenient numerical approach at the 
cost of requiring many more unknowns in the region external to the scatterer. 
An alternate technique is to couple the internal FEM solution to the exterior, 
unbounded region through an integral equation formulation provided by the 
boundary element method (BEM) [7]. The BEM provides an efficient solution 
to the wave equation in a homogeneous region. 

The FEM-BEM [8, 9, 10] method combines the strength of both of these 

methods. A single boundary is established at the surface of the material body. 

Inside the boundary, where the material properties may be inhomogeneous, 

the FEM solves for the fields. Outside the boundary, the BEM solution 

for free space is used which couples the incident field to the far field using 

the boundary conditions derived from the FEM. The link between the FEM 

and BEM is based on the boundary conditions of the differential equation 

at a discontinuity of the material properties. For electromagnetic fields, this 

— * 

reduces to continuity of tangential E and H across the surface of the material 
body. 

This section will show how the FEM and the BEM are coupled together 
to solve for both the scattered far fields and the internal fields for an incident 
plane wave on an inhomogeneous two-dimensional body. 

Differential Equation 

The scattering problem being considered here is two dimensional. The mate- 
rial properties of the body are constant along the z axis and the incident field 
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is also constant in z. Under these conditions, Maxwell’s equations separate 
into two equations which depend only on the values of the E,H fields in the 
z direction, 


and 



(2.31) 


(2.32) 


The relative permittivity £ r , permeability fi T , and wave number k = ^^//x r e r 
may depend on position x,y. The derivative V is in the x,y direction only. 
The other field components may be derived from E z and H x . 


For a transverse magnetic (TM) field, 


Ein c = E 0 ze~ jkiP (2.33) 

and the scattered field will also be TM. Likewise a transverse electric (TE) 
incident field, 


H„ 


H 0 ze 


Z a -jki-p 


(2.34) 


will produce only TE scattered fields. Both the TE and TM differential 
equations have the form of the scalar equation 

^V-(^V<£)+ =0 (2.35) 

with an incident field defined as 

<f>,nc = he-fc*. (2.36) 
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This scalar equation is the differential equation which will be solved by the 
FEM-BEM with the understanding that the solution corresponds to either 
the TM case with <f> = E 2 , v = //. r or to the TE case with <j> = H z ,v = e r . 


Boundary Conditions 

At a boundary between material 1 and material 2 with parameters v \ , k\ and 
V 2 ,k 2 respectively, the differential equation provides the following boundary 
conditions, 


4 > i — 4*2 


(2.37) 


1 d<f>j _ 1 d4> 2 
i/i dn v 2 dn 


(2.38) 


where h is an outward unit vector normal to the boundary and £ is a 
derivative in the n direction. The last boundary condition can be written 


V>i 


i>2 


l d4> 

* = ZK- 

The fields 4 > an d V’ are continuous across the boundary. 


(2.39) 

(2.40) 


Let t = z x n be a unit vector tangential to the surface. Then for TM 
fields 


4> = E 2 

f * — * 

4> = ju>t • H 

and for TE fields 

4> = H z 

xp = —jut • E. 


(2.41) 

(2.42) 


(2.43) 

(2.44) 
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The boundary conditions on <p and ip are equivalent to enforcing continuity 
of tangential E and H. 

FEM 

The FEM is ideal for solving the differential equation shown in Equa- 
tion (2.35). To construct a solution, a mesh of triangular elements is made 
which completely encloses the material body. As usual this mesh is made 
small enough to accurately model the inhomogeneity of the body and the 
variation of the fields. Nodes on the edges of the triangles represent the 
value of (p at that point. These nodal values are used to extrapolate a linear 
or quadratic solution in the element. Each node i gives rise to a generating 
function io, which extends across each element node t touches, has unit value 
at node i and has zero value at every other node. The total field can be 
written as a linear combination of these generating functions, 

N 

<P = 'Z<t>M ( 2 - 45 ) 

1=1 

where <p , is the value at node i for the N nodes on the mesh. 

The FEM solution seeks to minimize the residual of the differential equa- 
tion as weighted by each generating function Wj. The residual for the jth 
generating function is 

V • V<P ) + Wjds. (2.46) 

After an integration by parts, this residual becomes, 

Rj = — J — (Vwj • V<p — k 2 Wj<p) ds + J Wjipdl. (2-47) 
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The region of integration is S which is enclosed by the boundary C. The 
FEM boundary coincides with C. As before, ip = on the boundary and 
is assumed independent (known). There are N residual equations. Using the 
expansion in Equation (2.45) for <p and setting the residuals to zero gives N 
equations in the N unknowns 

r N r 1 

J uijipdl = ^2<f>i J - (Vu>j • Vt Vi — k 2 WjWi ) ds. (2-48) 

This is the matrix equation, 


Vj = Sj,<Pi (2.49) 

Sji = f — (Vt Vj ■ — k 2 wjWi) ds (2.50) 

Vj= [ WjiPdl (2.51) 

J G 

with Sji an N x N matrix and Vj an N dimensional vector. There is an 
implied summation over repeated indices in this and the following matrix 
equations. 

For a node k on the boundary, there is a generating function iu®, restricted 
to the boundary, which has unit value at node k and is zero at all other 
boundary nodes. These generating functions provide a means for expanding 
ip on the boundary, 

$ = '52'f>kW% (2.52) 

t=i 

in terms of ipk, the value of ip at the fcth boundary node. The total number 
of boundary nodes is Nb . Substituting this expansion for ip into the formula 
for Vj gives the matrix equation, 


Vj = Tj k ipk 


(2.53) 
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(2.54) 



where Tj k is a nonsquare matrix with N rows and Nb columns. The FEM 
matrix equation becomes, 


T ]k i> k = Sji<f>{. (2.55) 

Inverting the FEM matrix Sji and solving for fa gives, 

fa = (s-') tj T jk l> k . (2.56) 

This equation gives the internal field <f> in terms of ip, the normal derivative 
of <p on the boundary. Note that the index i runs from 1 to N for the total 
number of nodes while k runs from 1 to Ng for the boundary nodes. For the 
linkage to the BEM, cp need only be known on the boundary. We therefore 
restrict t to these boundary nodes with the index k' which runs from 1 to 

N b . 


4>k> = Z/t'/tV'fc (2.57) 

/« = (S-%^. (2.58) 

The matrix I k ' k is a square Nb x Nb matrix which relates <p on the bound- 
ary to xp on the boundary. This matrix acts like an impedance boundary 
condition and will be used as part of the BEM solution. 

In the formula for the impedance matrix, the inverse of the FEM matrix 
Sji is shown. Finding this inverse is not necessary. The matrix Sji is banded 
and large so that taking the inverse could be very costly. Given a value of 
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ip on the boundary, the solution of <p throughout the internal region and on 
the boundary can be obtained quickly using matrix algorithms for banded 
matrices. The mth column of the impedance matrix is therefore obtained by 
setting ip m = 1 and ipk = 0 for k ^ m and solving for <pk '• The mth column 
of the impedance matrix is then given by, 

Ik'm = <Pk '• (2.59) 

This method requires Ng iterations to build the Ng columns of . 

BEM 

The BEM is an integral equation method for solving fields in a homogeneous 
region. This integral equation is reduced using moment method techniques 
to a matrix equation which relates the known incident field to a set of un- 
known coefficients. Solution of the matrix equation solves for the unknown 
coefficients and from these coefficients, the scattered field may be determined. 

Consider a region S bounded by the curve C. This time, however, the 
area exterior to S will be included in the analysis. Because the exterior region 
is free space and homogeneous, v + = 1 and fc + = fcy = ^. The superscript -f 
refers to the exterior region. Inside S, the parameters v~ and k~ may vary. 

The integral equation is usually derived from Green’s second identity 
which yields the following relation for the field <p exterior to S in terms of 
the value of <p and ip on the boundary curve C, 

* = An, - f c Go(p,P)^{P) - (2.60) 

The outward normal is n' and ^ is the derivative along the normal. Since 
<p and rp are continuous across a material boundary, no superscript is needed 


20 



to distinguish the interior and exterior fields. The Green’s function solution 
for the two dimensional wave equation in free space is, 


g«(p,P) = W- Pll)- (2.6i) 

The zeroth order Hankel function of the second kind is consistent with the 
e> ul time convention. 

This integral relation is valid everywhere in the exterior region and can 
be used to find the far field once <p and ip are known. The integral equation 
results from taking the field point on the boundary. Rearranging terms, 

<P,nc = ^J c G l) (p,p')v> + 'lJj(p')dl' 

+ [* - l c SG °w pt) #?')■»' ■ (2-62) 

This equation relates the unknown <p and ip on the boundary in terms of the 
known incident field <pi„ c . Substituting the expansion of <p and ip on the basis 
of boundary weighting functions into this integral equation, 

»b r r 

<P.nc = HV’fc / G 0 (p,p')v + wj'(p , )dr 

k = i uc 

+ E h' [«? - J c ■ (2-63) 

The moment method testing functions are chosen to be the boundary weight- 
ing functions w? . This gives a Galerkin form of moment method solution 
since the expansion and testing functions are the same. Applying each of 
these testing functions to the integral equation gives the matrix equation, 

Vj = + Mf k Uk' (2.64) 

where 

Vj = f w?<p inc dl' (2.65) 

JC 
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( 2 . 66 ) 


m%* = j c J c wf(p)G4pJX(?)u’<a 


K t 


■=Uc W >" M ' il -Uc 


w 




8G„(p,?) 

dn' 


w%{p')dl'dl (2.67) 


The matrices Mj^ + , are TV's x TVb square matrices. Mj^ + is symmetric, 
but is not. There are Nb equations and 2 Nb unknowns in this matrix 
equation. Another set of Nb equations is needed to complete the solution. 


So far no information about the interior region has been used. If the 
interior is homogeneous, another integral equation can be derived which de- 
pends on the value of <p and ip on the boundary and uses the interior Green’s 
function solution with parameters v ~ , k~ . This provides another Nb set of 
equations in the 2 Nb unknowns, <pk> and ipk- Simultaneous solution of the 
resulting 2 Nb equations gives the BEM solution for the scattering from a ho- 
mogeneous dielectric body. The far field may be found using Equation (2.60) 
along with the asymptotic form of the Hankel function, 


<t>far(pp) = 


e -J*oP 





sff> 


X>^ + I w B k {p')e^'dV 

k = l Jc 


+ ^2 <Pk< I -jko{ri • r)w B (p')e ]k ° p p dl' 

k '= i Jc 


( 2 . 68 ) 


This equation gives the far field (p in the p direction in terms of the boundary 
fields at the nodal points <pk',ipk- 


For an inhomogeneous body, additional internal boundaries can be in-, 
eluded and corresponding integral equations derived. The number of un- 
knowns grows rapidly and can quickly become unmanageable. What is 
needed by the BEM is an additional relation between <p and on the bound- 
ary. This is exactly what the FEM interior solution provides. 
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FEM-BEM 


The FEM and BEM are combined by substituting the impedance matrix 
Equation (2.58) derived from the FEM into the BEM matrix Equation (2.64). 
This operation uses the continuity of (p and xp across the boundary curve C . 
The result is 


Vj = A jk xP k (2.69) 

A jk = Mf k + + Mjvh'k- (2.70) 

The unknowns in this matrix equation are the values of xp on the boundary. 
The solution for xp k , 

*‘ = ( A '% V ' < 2 - M > 

depends on the incident field through Vj. The matrix Aj k is independent of 
the incident field. 

Once the xp k are known, the internal field <p may be found by solving the 
FEM matrix Equation (2.56) with the known values of xp k . The impedance 
matrix is used to find the boundary value of <p as in Equation (2.58). The 
field can then be found in the far field by using Equation (2.68) and the <p, 
xp boundary values. 

Ill Matrix Inversion 

The resulting system of equations to be solved involves the inversion of a 
square matrix. This matrix is very sparse due to the limited coupling be- 
tween elements. The location of non-zero elements in the matrix is dependent 
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Figure 2.3: Matrix distribution of non-zero elements before renumbering. 


upon the nodal numbering. Numbering of the nodes should be done to al- 
ways minimize the difference of neighboring nodes. This allows the matrix 
to become very banded for greater inversion and internal memory storage 
efficiency. 


Presently, a renumbering algorithm [12] is used to renumber the- nodes. 
Figures 2.3 and 2.4 demonstrate the element distribution of a typical matrix 
composed with 123 nodes before and after the renumbering, respectively. The 
darkened squares represent non-zero elements. Note that the bandwidth of 
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Figure 2.4: Matrix distribution of non-zero elements after renumbering. 
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the matrix is dramatically reduced. The initial and final bandwidth was 70 
and 17, respectively for this particular matrix. The bandwidth is defined here 
to be the maximum difference minus one of the non-zero matrix elements for 
any row. 

The actual solver used employs simple Gaussian elimination to obtain the 
unknown nodal values. This appears to be adequate when the elements have 
approximately the same area size. When there is a drastic difference between 
element size, a partial pivoting algorithm would be required. However, this 
requirement would require more storage than presently used since only the 
banded elements of the matrix are stored. 


IV Scattered Far Fields 


Although the FEM solution provides fields in the vicinity of the scattering 
geometry, most often the scattered far field is of interest. A useful normal- 
ization quantity for the scattered far field is the echo width as defined below: 


Echo Width = 2tt lim 

p— * OO 


lltf'll * 

V‘ll 


(2.72) 


where U a represents the scattered field for either polarization. 

Obtaining the scattered far field for the ABC and BEM approaches is 
slightly different. The ABC approach uses the calculated total field at some 
specified circular boundary. The scattered field can be represented as 

V= £ aJ-HKKp)?"* (2.73) 


where o„ are coefficients and H„ are Hankel functions of the second kind. 
Once the coefficients are found, the scattered far field is obtained by using 
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the large argument form for the Hankel function. The scattered far field is 
then given by 


00 1 7 i p-jxp 

(274) 

The coefficients a n are determined by the following expression, indepen- 
dent upon how the total fields, U l , are calculated: 

3 n h /o 2w U'{a, 4)e**d4> - (-1)" -J n (Ka)e^ 

where U l is the total field at p = a and J„ is a Bessel function. The integra- 
tion in this equation yields the nth order contribution for the scattered and 
incident fields. The scattered field is then obtained by subtracting the nth 
order incident field (plane wave) contribution from the nth order total field. 

The BEM approach instead uses equivalent currents which are calculated 
at the coupling boundary surrounding the scatterer. The fields from these 
currents are then radiated into the far field using the standard radiation 
integral as discussed earlier. 

V Geometry Specification 

For any numerical solution, the geometry under consideration has to be spec- 
ified. The approach chosen here is to use a commercially available CAD 
package and mesher to provide the scattering code required data. Several 
packages are available at very reasonable costs. The package chosen here is 
provided by Structural Research and Analysis Corporation [4] and is called 
Geostar. They have implementations available on several computers ranging 
Rom PC’s to workstations. 


(2.75) 
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The scattering code reads an ASCII output file generated by Geostar 
which describes the geometry in terms of first or second order triangular 
elements. This file lists the coordinates of the nodes, the global node num- 
bering for each element, material specification for each element and required 
information to enforce the necessary boundary conditions. 
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Chapter 3 
Examples 


Two different example classes are presented. One is to validate the computer 
solution through comparisons with results generated by eigenfunction solu- 
tions. The other examines the scattered fields from some simple terminations. 
The calculated results shown here were all obtained using the FEM/BEM ap- 
proach. The FEM/ABC appears not as robust and requires a much greater 
mesh and hence more unkowns. The greater number of unknowns may not 
be as significant interms of compution time and required storage as the error 
resulting from the field propagating through the extra mesh. It is desirable 
to keep the meshed area as small as possible to minimize discretization error. 
Hence, the FEM/BEM is a natural for this. A future report will document 
the FEM/ABC performance. 

I Validation 

Code validation usually results through the comparison of its calculations 
to measurements or calculations through other independent solutions. An 
eigenfunction solution is considered to be the most exact form of numerical 
solution available. For this validation, comparisons will be made to two 
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Table 3.1: Geometry data 


Geometry 

Nodes 


Bandwidth 

Run Time (s) 


1577 

3016 

52 

126 

i 

90 

120 

9 

3.2 


2087 

3873 

45 

387 


2087 

3873 

45 

393 


1652 

3030 

41 

280 


geometries where eigenfunction solutions exist. The validation geometries 
will be (1) two concentric dielectric cylinders and (2) a perfectly conducting 
flat strip. The excitation in both cases will be from plane wave incidence as 
illustrated in Figure 3.1. 

The dimensions for the dielectric cylinder are 7 and 16 cm for the inner 
and outer radii, respectively. The relative permittivities for the inner and 
outer regions are e r = (4.,— j.2) and e r = (2.,—j.l). The relative perme- 
ability for both regions is unity. The width for the strip geometry is 10 cm. 
The frequency of the incident field for both cases is 3 GHz. The meshed 
geometries are illustrated in Figures 3.2 and 3.3. The mesh density used 
for both cases was approximately 10 nodes per linear length per free space 
wavelength. The BEM boundary was located at the perimeter of the meshed 
regions. Table I contains mesh information and computer run time for all 
the geometries. 

Figures 3.4 and 3.5 compare the TE and TM, respectively, FEM/BEM 
and eigenfunction bistatic echo width solutions for the concentric dielec- 
tric cylinders. Figures 3.6 and 3.7 compare the TE and TM, respectively, 
FEM/BEM and eigenfunction monostatic echo width solutions for the flat 
strip. All calculations used linear shape functions. Agreement between the 
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(a) 



p 10cm p 

(b) 

Figure 3.1: Test geometries: (a) Concentric cylinders, (b) Flat strip. 
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Figure 3.2: Mesh used for the concentric cylinders. 


32 



Figure 3.3: Mesh used for the flat strip. The centerline of the dashed region 
is the location of the strip. 

eigenfunction and FEM/BEM results improve as the mesh density is in- 
creased and/or higher order expansion functions are used. 

II Terminations 

Three examples are shown here to demonstrate simple triangular termina- 
tions. The terminations are lossy dielectric structures with a conducting 
backplane and are designated as Tl, T2 and T3. Terminations T1 and T2 
are shown in Figure 3.8a. and termination T3 is shown in Figure 3.8b. 
The relative dielectric constant for termination Tl had a uniform value of 
e r = (4.,— .05). Termination T2 had its relative dielectric constant linearly 
graded from c r = (2., —.05) at its tip to e r = (4., —.2) at its base. The relative 
dielectric constant for termination T3 had a uniform value of e r = (4., —.05). 
The relative permbeality for all three cases was fi = (l.,0.). 

Figures 3.9 through 3.11 contain the TM and TE echo widths for termi- 
nations Tl through T3, respectively. All calculations were performed at 3 
GHz with a linear nodal density of 20 per material wavelength. 
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Figure 3.5: TE bistatic echo width for the dielectric concentric cylinders. 
FEM/BEM - solid, Eigenfunction - dashed. 
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Figure 3.6: TM bistatic echo width for the flat strip. FEM/BEM - solid, 
Eigenfunction - dashed. • 


36 




Angle (Deg) 


Figure 3.7: TE bistatic echo width for the flat strip. FEM/BEM - solid' 
Eigenfunction - dashed. 
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(a) 



(b) 


Figure 3.8: Test terminations, (a) Terminations T1 and T2, (b) Termination 
T3. 
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Figure 3.9: Echo widths for termination Tl. TM - solid, TE - dashed. 
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Figure 3.10: Echo widths for termination T2. TM - solid, TE - dashed. 
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Figure 3.11: Echo widths for termination T3. TM - solid, TE - dashed. 
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A final plot is shown in Figure 3.12 which illustrates the internal TE 
scattered fields for termination T3 at an incidence angle of 25°. This angle 
was chosen due to the large echo width return, however, the correlation 
between near field distribution and echo width value is not simple. This 
plot does illustrate the standing wave pattern from waves due to localized 
scattering centers. These scattering centers can conceptually be identified by 
tracking the standing wave nature for a geometry when its electrical size is 
sufficiently large. The benefit of knowing these locations is that appropriate 
action can be taken to treat them from a scattering viewpoint. 
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Chapter 4 


Conclusions 


The status of a computer code for the calculation of the electromagnetic 
scattering for two dimensional geometries is presented. The code is based 
upon a finite element (FEM) solution for a bounded region. Two different 
approaches are examined to extend the solution for an unbounded region. 
One is using an absorbing boundary condition to simulate a reflectionless 
boundary (ABC) and the other one couples interior and exterior fields at a 
boundary through an integral equation (BEM). The advantages of the present 
FEM/BEM approach is that 1) no field approximations are made as in the 
ABC approach, 2) the boundary is placed at the outer contour of the geom- 
etry which minimizes the FEM unknowns and 3) this particular FEM/BEM 
approach, unlike others, retains the full advantage of the sparseness associ- 
ated with the FEM interior region for matrix inversion. 

Future activities include both the refinement of the existing code and 
the study of particular terminations of interest using both numerical and 
experimental approaches. 
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Appendix A 

BEM/FEM Numerical 
Implementation 

I Introduction 

This appendix will consider in more detail the calculation of the matrices 
Mjl + , Mf k t and the vector Vj used in the FEM-BEM. The relationship of Vj 
to the far field calculation will also be outlined. 


II Geometry 

The boundary for the FEM-BEM is defined by the mesh used in the FEM 
solution. This consists of a set of geometry nodes connected by line segments. 
The line segments form the boundary elements. Let the number of geometry 
nodes, which equals the number of boundary elements, be N g . Figure A.l 
shows a boundary with N g = 5. The boundary elements are indexed such 
that element number e is between geometry nodes e and (e mod N g )+ 1. The 
following discussion for the boundary element implementation includes both 
linear and quadratic generating functions. 
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For linear generating functions, the boundary nodes used as the unknowns 
for the FEM-BEM are the same as the geometry nodes, Nb = N g . Consider 
the set of linear generating functions for Figure A.l as shown in Figure A.2. 
The generating functions are combinations of the elemental linear functions 
JV* ■ (i = 1,2; e = 1 ,N g ) defined on each boundary element e as, 

*?.i(0 = ■!-* (A.l) 

*.%(0 = t (A.2) 

where £ is a normalized distance along the element (0 < £ < 1). For example, 
the linear generating function w 4 is given by, 


w 4 = Nl 2 + Nl v (A.3) 

For any linear generating function k, 

w h = A^( 2 ( fc _ 2 ) mo d7V ff +l),2 + Afjfc,!- (A.4) 

This formula gives the subscripts of Nl , as a function of k. The inverse 
relationship is also useful. Each function jV|, is a part of exactly one gener- 
ating function tu*. For linear elements, the kth boundary generating function 
which has JV^ ■ as part of its expansion is 

fc s= e (i = 1) 

k = (e mod ATg) -f 1 (t = 2) 


(A.5) 
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Figure A.3: Elemental linear generating functions on a single boundary ele- 
ment. 


The set of quadratic generating functions for Figure A.l is shown in Fig- 
ure A.4. The interior boundary nodes are shown as open circles. The numbers 
correspond to the boundary elements. These generating functions are com- 
binations of the elemental quadratic functions iV,?,- (i = l,3;e = 1, N g ). The 
three functions defined on boundary element e are, 


KdO = 2f * -H + 1 (A.«) 

W 1 W 2 W 3 W 4 W 5 W 6 W 7 W 8 W 9 W !0 W 1 



Figure A.4: Quadratic generating functions on the boundary elements. 
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0 1 


Figure A. 5: Elemental quadratic generating functions on a single boundary 
element. 

KM) = ~4£ 2 + 4£ (A.7) 

^(0 = (A.8) 

where £ is again the normalized distance along the boundary element (0 < 
( < 1). The expansion of wj on the set of functions JV 3 ,- depends on whether 
j is even (an internal node) or odd (a geometry node), 

W 2jl = Nj, 2 (A.9) 

w 2j'-l = ^ ( 3 (j -,_ 2)modA r ff+1) 3 + Nj, a (A. 10) 

for j' = l,N g . Notice that W 2 y extends over one boundary element while 
w 2 j>-i extends over two. Each function JV 3 ,• is a part of exactly one quadratic 
generating function. The kth quadratic boundary function which has JV 3 , as 
part of its expansion is 

k = 2e— 1 (i = 1) 

k = 2e (t = 2) (A. 11) 

k = 2(e mod N g ) + 1 (t = 3). 
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The matrices Mj^ + and Mj k t are defined by, 


II 

V 

JcJc p') w ? (p') dl ' dl 

(A. 12) 

j§; 

ii 

j c j c J c ™?{p) ^w*,(P)dl dl. 

(A.13) 


The actual integrations will be done over the functions N^ t and Nj n resulting 
in the terms 


- " + } f NUi)G^?)KAC)it' 


(A. 14) 


m*, = Sc N tM) N !AO dl ~ J c j c KAO 8G " a ( * - KAO<*dH A-is) 


These integrals are defined over straight line elements. For a given pair of 
•boundary elements e and /, the integrals for (i = l,p) and (j = l,p) can be 
done simultaneously thus saving run time in function calls for the Green’s 
function. The terms m^. e i and my+ ;e , are added to the appropriate matrix 
elements of M fp and based on the formulas in Equations (A.5) or 

(A. 11). 

The voltage vector Vj is defined as 

Vj = j c w*<f> tnc dl'. (A. 16) 

The integrals 

V,,i = J c KAO*incdl' (A. 17) 


are calculated for each boundary element e with (t = l,N g ) and combined 
to form the Vj vector. 
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The following notation and assumptions will be made for the BEM in- 
tegrations. The integral around the boundary will be done in a counter- 
clockwise direction. This defines the integration path along each boundary 
element. Let p e be the position of the geometry node of boundary element 
e where the integral begins. The vector from one end point to the other is 
l e where l e =|| l e || is the length of boundary element e. Any point on the 
boundary element is given by 

P = P* + & (A. 18) 

for 0 <*<1. The outward normal and the exterior region will always be to 
the right of the integration path. The outward unit normal is then given by 

ti e = (l e x zj /l e . (A. 19) 


III Component 

The free space Green’s function for the two dimensional wave equation is, 

a„[p,P) = f\k II f - ? II). (A.20) 

Using this Green’s function in Equation (A. 14) for the term rnj+. ei and 
parameterizing the line integrals with Equation (A. 18) gives, 

m Z,.i = /' /' N UOKM')B^\k || p- ? \\)U'H (A. 21) 

P = P! + U (A.22) 

? = P,+ W (A.23) 
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e - P = (p> - ft) + ('"/{ - «')• 


(A.24) 


From the symmetry of this integral, 


m 


<J>+ 

e,«;/,n 


= m 


/.n;e,i" 


(A. 25) 


This means that we only need to do a calculation of m^ ; ■ for e < /. The 
symmetry gives all the cases for e > /. 

The Hankel function of order zero has a logarithmic singularity at || 
p — p' ||= 0. When elements e and / are separated, the integrand is well 
behaved and a simple Gaussian quadrature is sufficient to calculate the dou- 
ble integral. When elements e and / touch at an endpoint, there is one 
logarithmic singular point in the unit square of £,£' space where the inte- 
gration is done. This is a very mild singularity and can be ignored when a 
Gaussian quadrature is used. 

When e = /, the logarithmic singularity becomes a line (£ = £') in the 
unit square. This double integral must be done more carefully to get accurate 
answers. For e = /, 

= z ~ ± £ KnU) {£ if- cm') <jf.(A. 28 > 

The innermost integral will be examined. It has the form, 

{> = fnm' 

Jo 

= £ new + £ new 


(A.27) 
(A. 28) 
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where /(£') has a log singularity at The singularity will be moved to 

zero in the integration variable with the transformation t = in the first 
integral and t = in the second, 

{ } = fa m - o) + (i - 0 / (* + (i - oo dt. (a.29) 

J u 

Substituting the integrand of the double integral for /(£'), 

{} = iff 1 (*«0 

J o 

+(1 - i)K< (( + (1 - 00 (*1,(1 - 00 *■ (A.30) 

Let N * ,(*) = a 2 a: 2 + a i* + Go- This equation works for linear ( p = 2) or 
quadratic (p = 3) functions. Then, 

Nc,i (£(1 — 0) = ^ a 2^ 2 + t (~ — 2 a 2 £ 2 ) + (a 0 + a.i£ + a 2 £ 2 ) (A. 31) 

^e,i (£ + (!“ 00 = ^ 2 ® 2(1 - 0 2 + * ( a l(l - 0 + 2a 2 £(l - 0) 

~h( a u + + fl20)- (A. 32) 

Substituting this expansion in powers of the integration variable t into Equa- 
tion (A.30) gives, 

{> =( 

-f(a 0 + aj£ + a 2 0)^o(^fcO} 

+(1-0 {a 2 (i-0 2 MM e (i-0) 

+ ( ai (l - 0 + 2a 2 ^(l - 0) h (fc/e(l - 0) 

+ (a 0 + a,\t + a 2 0)^o (fcl e (l — 0)} (A.33) 


where the /*(<*) are defined as 



The reduction of these integrals was done in [8, 11] and is repeated here. 

Not much can be done with / 0 (a) -except integrate the series expansion 
of 






(A. 35) 
(A. 36) 


- | (in (§) + 7 ) M,) (jft i (A.37) 

where 7 = .57721 ... is Euler’s constant. Grinding through the integral gives 

k 


[ Z H {2) (z'\dz' = ? V' ^ ( —\ * 

'» u ( 1 £(«)* Uj 2*+i 

.[ i _ J '2{ in (i ) + 7 __L T -|;i} 

/«(«)= - r H^\z')dz'. 

a Jo 


(A. 38) 


(A. 39) 


The last summation is defined to be zero when k = 0, X)«=i - = 0. The 
integral for /i(a) can be done exactly, 


h(z\a) = j 7]H^ 2) (arj)drj 

= is:^^ 


(A.40) 
(A.41) 
(A. 42) 


Using the small argument approximation of iff '(0:7 ) at rj cs 0 gives, 


/,(«) = - 
a 


1 r 




(A.43) 
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/ 2 (a) is found by an integration by parts, 


I 2 {z\a) = [ rj 2 H^\a , q)d , q 

J 0 

(A.44) 

= / Vl[{v\<*)dv 

Jo 

(A.45) 

= vhiv;*) lo - / 

(A.46) 

= vli (»?; <*) lo + ^ J o (tfo 2) (a*?)) dr] 

(A.47) 

= Vliivi*) IS +A 1 o - f Ho\av)dv 

a' l yo J 

(A.48) 

/ 2 (a) = 1//PV) + AfffV) - -M a )- 

a or a' 

(A.49) 


The e = / term is evaluated by doing a Gaussian quadrature on the 
outer integral in Equation (A. 26). For each value of f, the inner integral is 
evaluated using Equation (A. 33) and the /jt(a) functions given above. 


IV Component 


The normal derivative of the free space Green’s function is, 

dG „(p,P) 


dn' 


= n' • VG 0 {p,p') 

= z£ H '?\k\\p-t\\)n'.R 


(A. 50) 


where R is a unit vector from the source to the field point, 


R = (p - ?)/ II p- P' II • 


(A.51) 


When e = /, there is a singularity in the second integral in Equation (A. 15) 
for mt This singularity comes from the ; behavior of the Hankel 
function of order one. Also, with e = / on a linear path of integration, n'’R = 
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0 except at p = p'. The net effect of the second integral in Equation (A. 15) 


is to reduce the first integral by one half resulting in, 

"£n = \j c KM)KAO<u 

= \ £ N UQKMW- (A.52) 

This integral can be done exactly which gives 

[ KM)KAm = £ KAWUM = \ (a.s3) 

£ KMWUm = J (A.54) 

for linear generating functions and 

1‘ KAOKMW = /' KAOKAM = ^ (a.m) 

/' NUWUtW = Jg (A. 56) 

£ *£,(()*&«)<« = / o ' NlMftlAM = ~ (A.57) 

’ jf KlUWMW = ^ (A.58) 

for quadratic generating functions. 


When e//, Equation (A. 15) reduces to, 


„,<H- 


jkljj 


f f II a-p' ll)(< • 

Jo Jo 



(A. 59) 
(A.60) 


P’ = Pe + fe{' 

p-p> = (p> - /0 e ) + (l/£ _ l e ^')- 


(A.61) 
(A. 62) 


In general »n££. e>f - ^ m tXf,n because of the term (to' • i2). However, by cal- 
culating both (to' • R ) and (n'^ • J?), both terms - and *nf+.j tn can be 

calculated in the same integration loop. 


56 



For e/ / and e not touching /, the integrand of Equation (A.59) is well 
behaved and can be found by a simple Gaussian quadrature. 


When e ^ / and e touches / at a corner, there is a singularity at that 
corner. While this singularity is not severe, it can be easily removed. We 


assume f — e mod N g + 1 so that pj = p e + l e . Let 17 = 1 — (' ( d-rj = — d£') 
in Equation (A.59), 

/„' SI nuonza 1 - ») 

II A - P ll)(< ■ (A.63) 

- P = P,+ f/< (A.64) 

p' = pj -l e i) (A. 65) 

p-p' = (A. 66) 


The double integral over the unit square in £,77 space is divided into two 
triangular regions as shown in Figure A. 6, 


TO 


d>+ 

f,n;e 


= TO 


upper 

/,n;e,i 


-(- TO 


lower 


(A.67) 


For the lower triangle, let £ = t, rj = th. This gives 

"fc = ^ /„’ IS N U‘)K •(! - <*) ' 

|| p-p' ||)(ra' • (A. 68) 

p-p' = t{l f + l e h). (A. 69) 


Notice how the singularity at t = 0 has been removed since the term tH[ 2 \k || 
p — p' ||) is bounded at t = 0. For the upper triangle, let rj = s, { = sw. This 
gives 


TO 


upper 

/,n;e,t 


j| p — p' ||)(n' • R)sdsdw 


(A.70) 



V 


w 



0 t I f 


Figure A. 6: The division of the unit square into upper and lower triangles. 
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The singularity at s = 0 has been removed. Since the integrands for tti'j’H 
and are well behaved, a Gaussian quadrature is used for the double 

integrals. Adding the results as in Equation (A. 67) gives m^. e i for the case 
e < / and touching. The case e > f and touching is found by calculating 
the two terms m^ ;e , and simultaneously as indicated earlier. 

V Component v e j 

The induced voltage depends on the incident field. We assume here that 
the incident field is a plane wave of unit magnitude coming from the pj nc 
direction, 


Jy. p -jkinc‘P 

Y 'tnc — e 

(A.T2) 

■ — gjkoPiric'P 

(A. 73) 

Equation (A. 17) gives the term v e ,,, 

u e , ( = i e f 

Jo 

(A. 74) 

= l^jkoPinc fi' [' Nltt)e M d( 
Jo 

(A. 75) 

(L = fcoPinc ' 

(A.76) 

Substituting the polynomial expansion Nl(() 

— a2 ^ 2 4 - a.i( + do into this 

equation gives, 

v e< i — A e (02^2(0) + otJi(a) + aoJ(j(a)) 

(A. 77) 

A*, _ l^jkoPinc Pe 

(A.78) 

Ma) = /' 
Jo 

(A.79) 
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The Jfc(a) integrals can be reduced, 

Ma) = ±(e*-l) 

ja v ' 

Ji(a) = (e ja - Ma)) 
ja ' ' 

Ma) = ±-(e ja -2J 1 (a)) 
ja ' 7 


(A. 80) 
( A.81) 
(A. 82) 


for o / 0 . For a close to or equal to zero, the Jfc(a) integrals are best 
evaluated using the series expansion for e ja ^ , 


T (n.\ = V' 0 a ) 


(A. 83) 


VI Far Field 


The far field <f>f ar in the p/ ar direction in terms of the boundary fields tj>k > , iftk 
at the nodal points is, 

~n b 


p-}kof> ( ~ I O,’ \ n b r _ 

= —F-(j\Hr) £ / w?(p>y t °>‘->'di' 

\P \ TTKo / 'C 


+ 4>u f ~3^( n ' • Pf°r)wE(p')e 

k '= i *' c 


JkoP/arP (ft' 


.(AM) 


This can be written as, 


<f>far(pf„rP) = C (V’fc-F* + <f>k'F*S) 

Ft = [ w^(p')e iko ^'dl' 

* c 

* c 


(A.85) 

(A. 86) 
(A.87) 



(A. 88) 
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where an implicit summation over repeated indices has been used. Here 


again, the terms 


ft, = l /' 

J 0 

(A. 89) 

ft, = -iM"' • Ml [ 

(A. 90) 


are calculated for each boundary element e (t = l,p) and combined to give 
Fjf, F*. The factor (n' • pj ar ) is constant over the linear boundary element 
e and so it was removed from the integrand. 

By comparing the far field Equations (A. 89) and (A. 90) with the voltage 
Equation (A. 74), the far field terms Ff ti , can be written as, 


Fgj — A' c (a 2 J2(ffl ) + &i«fi(<0 + duJo(a')) (A. 91) 

F*i = -jko(n r ■ Pf ar )A' e (a 2 J 2 (a') + oj Ji(a') + a 0 Jo(a')) (A.92) 

where A' e and a' are defined by 

A' e = l e e jk (A. 93) 

a = k 0 p far • l e . (A.94) 


The J*. (a) functions were defined earlier. 


61 



Bibliography 


[1] M. Kragalott and E. Newman, “A User’s Manual for the General Cylin- 
der Code (GCYL),” Technical Report 722644-1, ElectroScience Labora- 
tory, The Ohio State University, June 1990; prepared for Rosemount, 

Inc, Burnsville, MN, under Contract No. PO D55951. 

[2] P. Silvester and R. Ferrari, Finite Element for Electrical Engineers, 2nd 
Edition, 1990 Cambridge University Press. 

[3] B. McDonald and A. Wexler, “Finite-Element Solution of Unbounded 
Field Problems,” IEEE Trans. Microwave Theory Tech., Vol. MTT-20, 
pp. 841-847, Dec. 1972. 

[4] Structural Research and Analysis Corporation, 1661 Lincoln Blvd, Suite 
200 Santa Monica, CA, 90404. 

[5] A. Bayliss and E. Turkel, “Radiation Boundary Conditions for Wavelike 
Equations,” Comm. Pure Appl. Math, Vol. 33, pp. 707-725, 1980. 

[6] D.J. Richards, A. Wexler, “Finite-Element Solutions within Curved 
Boundaries,” IEEE Transactions on Microwave Theory and Techniques, 

Vol. MTT-20, No. 10, October 1972. 

[7] J, Yashiro and S. Ohkawa, “Boundary Element Method for Electromag- 
netic Scattering from Cylinders,” IEEE Antennas Prop., Vol. AP-33, 

No. 4, pp. 383-389, April 1985. 

[8] A. Cangellaris and R. Lee, “The Bimoment Method for Two Dimensional 
Electromagnetic Scattering,” IEEE Antennas Prop., Vol. AP-38, No. 9, 
pp. 1429-1436, Sept. 1990. 

[9] K. Wu, G. Delisie, D. Fang, M. LeCours, “Coupled Finite Ele- 
ment and Boundary Element Methods in Electromagnetic Scattering,”, 

Finite Element and Finite Difference Methods in Electromagnetic Scattering , 
Ed. M. Morgan, Elsevier Science Publishing Co., 1990. 


62 




[10] J. Jin and V. Liepa, “Application of Hybrid Finite Element Method to 
Electromagnetic Scattering From Coated Cylinders,” IEEE Antennas 
Prop., Vol. AP-36, No. 1, pp. 50-54, Jan. 1988. 

[11] M. Koshiba and M. Suzuki, “Application of the Boundary-Element 
Method to Waveguide Discontinuities,” IEEE Microwave Theory Tech., 
Vol. MTT-34, No. 2, pp. 301-307, 1986. 

[12] R. Collins, “Bandwidth Reduction by Automatic Renumbering,” In- 
ternational Journal for Numerical Methods in Engineering, Vol. 6, pp. 
345-356, 1973. 


63 



723224-3, “Analysis of Lossy Composite Terminating Structures” Andre, Dominek, Munk and Wang, September 1991 



